/*  scfcatan.cpp  by K.Tsuru  */
// function ID = 9116  since ver 2.18
/****************************************************
SComplex class
It returns arctan(z). z = x+iy
Catan(z) = (i/2) * Clog{(i + z)/(i - z)}.
Catan(iy) = i*Atanh(y)

Let z = x [+|-] i*y
arctan(z) = arctan(x [+|-] i*y) ([+|-]y >= 0)
  1            2x              i             2y
=--- arctan ----------- [+|-] ---artanh-----------
  2          1-x^2 -y^2        2         1+x^2+y^2
*****************************************************/
#ifndef SN_H
#include "sn.h"
#endif
#define CatanIZMethod 1 // 25.0sec ...fastest
#define CatanXYMethod 0 // 25.9sec
   // WolframMethod     // 36.9 
static const SComplex I_2(0.0, 0.5); // i/2

SComplex Catan(const SComplex& z){
    if(z == 0.0) return 0.0;
    if(Re(z) == 0.0){ // pure imaginary z = iy
        if(Dabs(Im(z)) >= 1.0) z.Real().SetError(z.Real().DOMAIN_ERR, "Catan z = iy(|y| > 1.0))", 9116);
        return SComplex( 0.0, Atanh(Im(z)) ); // i*atanh(y)
    }
    SComplex r;

#if CatanIZMethod // faster than below
    SComplex p = (IU + z)/(IU - z);
    r = I_2 * Clog(p);

#elif CatanXYMethod
    SDouble x(Re(z)), y(Im(z)), zsq = x*x+y*y, d1 = zsq - 1.0;
    if(d1.IsHidden()) r = SComplex(MPi4(), Atanh(Sin(x))/2); // |z|=1
    else {
        SDouble p, q, d;
        p = 2*x; d = 1-zsq; q = Atan2(p, d); // instead of Atan(p/d) you must use Atan2(x, y)
        x = DsDiv(q, 2);

        p = 2*y/(1+zsq); q = Atanh(p); y = DsDiv(q, 2);
        r = SComplex(x, y);
    }
#else
// WolframMethod
// arctan z = (i/2)[ln(1-iz)-ln(1+iz)]
// 1-iz = 1-i(x+iy) = 1+y -ix, 1+iz = 1-y+ix
    SDouble x = Re(z), y = Im(z);
    SComplex p(1+y, -x), q(1-y,x); 
    r = Clog(p) - Clog(q);
    r = I_2 * r;
#endif

    return r;
}